#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <algorithm>

#include "Vector/vector.03.h"
#include "Color/color.01.h"
#include "Timer/timer.01.h"
#include <vector>
#include <deque> 

#define SQR(x) ((x)*(x))

#include "Data.h"


using namespace std;

class Kernel {
public:

	typedef unsigned int uint;
	typedef long unsigned int luint;
	typedef double real;
	
	Data &data;
	
	Kernel(Data &data) : data(data),
	relaxationF(false),
	potentialF(false),
	Vcl(0.99)
	{
		
		//pow13i7_i3 = pow(13.0/7, 1.0/3);
		
	}
	
	
	
	// grid
	class Grid{
		public:
			
			int xn, yn, zn;
			double xs, ys, zs;
		
			vector< vector<Data::P*> > node;
			
			inline vector<Data::P*> &operator()(int xi, int yi, int zi) {
				return node[xi + xn*(yi + yn*zi)];
			}
	
	} grid;
	
	// integrator
	double cdt;
	double dtmax;
	double dRmax;
	double dEmax; // minimum potential well depth

	double dt;
	double period;
	double ss_period;
		
	double g; 
	
	bool relaxationF;
	//bool gravityRelaxationF;
	bool potentialF;
	double Frx;
	double Vcl;
	
	//double pow13i7_i3;
	static const double pi05 = .5*3.141592654;
	
	string outputfolder;
	inline void calc_force(vector<Data::P*> &pi) {
	
	
		
		for (int i=0;i<pi.size()-1;i++)
			for (int j=i+1;j<pi.size();j++) {
						
					real dx = pi[i]->R.x - pi[j]->R.x; // 1
					real dy = pi[i]->R.y - pi[j]->R.y; // 1
					real dz = pi[i]->R.z - pi[j]->R.z; // 1
								
					real r2 = dx*dx + dy*dy + dz*dz; // 5

					Data::PW &pwij = data.pw[pi[i]->type][pi[j]->type];
					
					//real rc2 = pwij.rc2;
					
					if ( r2 > pwij.rc22 ) continue; 
					
					//real eij = data.pw[pi[i]->type][pi[j]->type].e;
					//real sij2 = data.pw[pi[i]->type][pi[j]->type].s;
					
					//real r12 = data.pw[pi[i]->type][pi[j]->type].rfdot02;//pow13i7_i3*sij2;
						
					//real o = (r12 < r2) ? ( sin(M_PI*((r2-r12)/(rc2-r12)+0.5))+1)*0.5 : 1.0;
					
					//real sij4 = sij2*sij2;
					//real sij6 = sij2*sij4;
					
					real r4 = r2*r2; // 1
					real ir6 = pwij.s6/(r2*r4); // 2
					real ir12 = ir6*ir6;  // 1

					real o = pwij.rc12 < r2 ? (sinf((r2-pwij.rc12)*pwij.pirc22rc12 + pi05)+1)*0.5 : 1.0;					
					real f = pwij.e12*ir6*r4*(ir12-ir6)*o; // 4
					
							
					real dxf = dx*f;
					real dyf = dy*f;
					real dzf = dz*f;
							
					pi[i]->F.x += dxf; // 2
					pi[i]->F.y += dyf; // 2
					pi[i]->F.z += dzf; // 2
						
					pi[j]->F.x -= dxf; // 2
					pi[j]->F.y -= dyf; // 2
					pi[j]->F.z -= dzf; // 2
					
					
					if (!potentialF) continue;
					
					real pot = pwij.e*(ir12 - .5*ir6);
					
					pi[i]->pot+=pot;
					pi[j]->pot+=pot;
					
							
			}
			
	}


	inline void calc_force(vector<Data::P*> &pi,  vector<Data::P*> &pj) {
					
		for (int i=0;i<pi.size();i++)
			for (int j=0;j<pj.size();j++) {
						
					real dx = pi[i]->R.x - pj[j]->R.x; // 1
					real dy = pi[i]->R.y - pj[j]->R.y; // 1
					real dz = pi[i]->R.z - pj[j]->R.z; // 1
								
					real r2 = dx*dx + dy*dy + dz*dz; // 5

					Data::PW &pwij = data.pw[pi[i]->type][pj[j]->type];
					
					if ( r2 > pwij.rc22 ) continue; 
					
					
					real r4 = r2*r2; // 1
					real ir6 = pwij.s6/(r2*r4); // 2
					real ir12 = ir6*ir6; 
					

					real o = pwij.rc12 < r2 ? (sinf((r2-pwij.rc12)*pwij.pirc22rc12 + pi05)+1)*0.5 : 1.0;
					real f = pwij.e12*r4*ir6*(ir12-ir6)*o; // 4
			
					real dxf = dx*f;
					real dyf = dy*f;
					real dzf = dz*f;
			
					pi[i]->F.x += dxf; // 2
					pi[i]->F.y += dyf; // 2
					pi[i]->F.z += dzf; // 2
			
					pj[j]->F.x -= dxf; // 2
					pj[j]->F.y -= dyf; // 2
					pj[j]->F.z -= dzf; // 2
					
					
					if (!potentialF) continue;
					
					real pot = pwij.e*(ir12 - .5*ir6);
					
					pi[i]->pot+=pot;
					pj[j]->pot+=pot;
					
				}
					
		}

	double min(double value1, double value2, double value3, double value4) {
			return std::min(std::min(std::min(value1, value2), value3), value4);
	}

	
	virtual void run(){
			
		Data &d = data;
		vector<Data::P> &p = data.p;
		
		for(int i=0; i<data.p_size; i++) d.p[i].F*=0;

		double t(0);
		double dt = 1e-5/d.tU; // 

		g = -9.81*SQR(d.tU)/d.lU;
	
		double Vmax(0);
		double Fmax(0);
	
		double Sf20; 
		
		int step(0);
	
		save_ss(-1);
		
		double ss_t = ss_period/d.tU;
		
		Timer tr;
		Timer trg;
		Timer trf;
		
		
		
		double grid_05lenx = 0.5*grid.xn*grid.xs;
		double grid_05leny = 0.5*grid.yn*grid.ys;
		double grid_lenz = grid.zn*grid.zs;
		
		grid.node.resize(grid.xn*grid.yn*grid.zn);
		
		
		do {
		
			tr.set();
			
			double dt05 = 0.5*dt;
				
			trg.set();
				

			for (int i=0; i<grid.node.size(); i++) 
				grid.node[i].resize(0);
				
			//piano_hilbert_ordering();

				
					
			for(int i=0; i<data.p_size; i++) {
			
					double dt05_over_m = dt05/d.pt[p[i].type].m;
					p[i].R+= (p[i].V +  dt05_over_m*p[i].F)*dt;
					p[i].V+= dt05_over_m*p[i].F;
				
					// boundaries	
					if (p[i].R.z < 0) {
						p[i].R.z*=-1;  p[i].V.z*=-1; 
					} 
					
					if (grid_05lenx < fabs(p[i].R.x) ||  grid_05leny < fabs(p[i].R.y) || grid_lenz < p[i].R.z) {
						p[i--] = p[--d.p_size];
						continue;
					} 
		
					// greed update
					grid(int((p[i].R.x+grid_05lenx)/grid.xs), int((p[i].R.y+grid_05leny)/grid.ys), int(p[i].R.z/grid.zs)).push_back(&p[i]);
					
					// forces
					p[i].F.x=p[i].F.y=p[i].F.z=0;
				}
				
				trg.get();

				trf.set();
				
				calc_all_forces();
				
				trf.get();

				double Sf2(0);
				double V2max(0), F2max(0);
				for(int i=0; i<d.p_size; i++) {
			
					
					if (!relaxationF) d.p[i].F.z += g*d.pt[p[i].type].m;
					
					d.p[i].V+= d.p[i].F*(dt05/d.pt[d.p[i].type].m);
					
					double f2i = p[i].F.sqr();
					double v2i = p[i].V.sqr();
										
					if (v2i > V2max) V2max = d.p[i].V.sqr();
					if (f2i > F2max) F2max = d.p[i].F.sqr();
					
					Sf2+= f2i;

				}
			
				t+=dt;
			
				
				if (relaxationF) {
					t=0; 
					for(int i=0; i<d.p_size; i++) p[i].V*=Vcl;
					if (step==0) Sf20=Sf2;
					if (Sf2/Sf20 < Frx) {
						relaxationF = false;
						remove_debries();
						save_ss(t);
					}
				}

				
				tr.get();		
		
				if ( fabs(t - ss_t) < 1e-12 ) { 
					save_ss(t*d.tU);
					ss_t+=ss_period;
					dt = dtmax;
				}
		
				step++;
		
				Vmax = sqrt(V2max);
				Fmax = sqrt(F2max);
		
				dt = min(dRmax/Vmax, dEmax/(Fmax*Vmax), cdt*dt, dtmax);
				//dt = 0.00001;
				
				if ( ss_t < t+dt) dt = ss_t-t;
		
				//if (step<300)
				
				//fprintf(stdout,"%9i  %10.5e    %8.2e    %8.2e  %8.2e    %6.2f  %6.2f   %6.2f    %10.2e, %10.2f", 
				//		  step, t,  dt, dRmax/Vmax, dEmax, trg.sin, trf.sin, tr.sin, d.p_size/tr.sin, tr.sin_tot);
				fprintf(stdout,"%9i  %10.5e   %7.2e   %8.2e  %8.2e    %10.2f  %10.2e, %10.2f",  step, t, Sf2/Sf20, dRmax/Vmax, dEmax/(Fmax*Vmax), tr.sin,  d.p_size/tr.sin, tr.sin_tot);		
				//fprintf(stdout,"%9i  %10.5e    %8.2e  %8.2e    %10.2f  %10.2e, %10.2f",  step, t, dRmax/Vmax, dEmax/(Fmax*Vmax), tr.sin,  d.p_size/tr.sin, tr.sin_tot);
				fflush(stdout);
				fprintf(stdout,"\n");
					
		} while (t<period);
	
		
	
	}

	
	inline void calc_all_forces() {
			
				for (int zi=0; zi<grid.zn; zi++)
				for (int yi=0; yi<grid.yn; yi++)
				for (int xi=0; xi<grid.xn; xi++) {
							
					vector<Data::P*> &pi = grid(xi, yi, zi);
			
					if (pi.size()==0) continue;
							
					calc_force(pi);
			
					
					if (0 < xi && 0 < yi && zi+1<grid.zn)					calc_force(pi, grid(xi-1,yi-1,zi+1)); // 12
					if (0 < xi && zi+1<grid.zn)							   calc_force(pi, grid(xi-1,yi  ,zi+1)); // 1			
					if (0 < xi && yi+1<grid.yn)								calc_force(pi, grid(xi-1,yi+1,zi  )); // 13
					if (0 < xi && yi+1<grid.yn && zi+1<grid.zn)			calc_force(pi, grid(xi-1,yi+1,zi+1)); // 6
			
					if (0 < yi && zi+1<grid.zn)								calc_force(pi, grid(xi  ,yi-1,zi+1)); // 4
					if (zi+1<grid.zn)												calc_force(pi, grid(xi  ,yi  ,zi+1)); // 2
					if (yi+1<grid.yn && zi+1<grid.zn)						calc_force(pi, grid(xi  ,yi+1,zi+1)); // 5
					if (yi+1<grid.yn)												calc_force(pi, grid(xi  ,yi+1,zi  )); // 8
			
					if (xi+1<grid.xn && 0 < yi && zi+1<grid.zn)			calc_force(pi, grid(xi+1,yi-1,zi+1)); // 11
					if (xi+1<grid.xn && zi+1<grid.zn)						calc_force(pi, grid(xi+1,yi  ,zi+1)); // 3
					if (xi+1<grid.xn )											calc_force(pi, grid(xi+1,yi  ,zi  )); // 10
					if (xi+1<grid.xn && yi+1<grid.yn && zi+1<grid.zn)	calc_force(pi, grid(xi+1,yi+1,zi+1)); // 7
					if (xi+1<grid.xn && yi+1<grid.yn)						calc_force(pi, grid(xi+1,yi+1,zi  )); // 9
					
				}
	}
	

	void remove_debries() {		
		
		printf("start to search and remove debries..."); fflush(stdout);

		for(int i=0; i<data.p_size; i++) 
			data.p[i].pot=0;
		
		potentialF = true;
		calc_all_forces();
		potentialF = false;
		
		for (int i=0; i<data.p_size; i++)
			if ( fabs(data.p[i].pot) < 10000 )
				data.p[i--] = data.p[--data.p_size];	
		
		printf("done\n");
	}
	
	
		
	virtual void save_ss(double t) {
	
		mkdir(outputfolder.c_str(), 0777);
			
		char filename[256];
		sprintf(filename, "%s/ss.%.2f.dat", outputfolder.c_str(), t);
		
		fprintf(stderr,"%s\n", filename);
			
		FILE *fp = fopen (filename , "wb");
		
		struct MEM { Vector<float> R; float type; } *mem;
		mem = new MEM[data.p_size];
	
		for (int i=0; i<data.p_size; i++) {
			mem[i].R = data.p[i].R;
			mem[i].type = data.p[i].type;
		}
		
		fwrite(mem, sizeof(MEM), data.p_size, fp);
		fclose(fp);
		delete [] mem;
	
	}
		
};


